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We study the implications of adopting hyperbolic driver coordinate conditions motivated by ge- 
ometrical considerations. In particular, conditions that minimize the rate of change of the metric 
variables. We analyze the properties of the resulting system of equations and their effect when 
implementing excision techniques. We find that commonly used coordinate conditions lead to a 
characteristic structure at the excision surface where some modes are not of outflow-type with re- 
spect to any excision boundary chosen inside the horizon. Thus, boundary conditions are required 
for these modes. Unfortunately, the specification of these conditions is a delicate issue as the out- 

- - - flow modes involve both gauge and main variables. As an alternative to these driver equations, 

we examine conditions derived from extremizing a scalar constructed from Killing's equation and 

^5 , present specific numerical examples. 
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The choice of suitable coordinate conditions has certainly played a major role in the understanding of solutions 
^ of Einstein equations at the analytical level. From the early well posedness result by Choquet-Bruhat 0| to recent 

■ global existence proofs a judicious choice of coordinates was key to yielding a tractable problem. 

I ' At the numerical level, coordinates play even a more crucial role as an unfortunate choice might render the simulation 
^ ■ useless despite much computational effort. In fact, the optimal situation is one where coordinates not only behave well 
fS| (i.e. not forming coordinate singularities) but also aid in the simulation. The latter refers to a choice of coordinates 
On . that adapts to the problem at hand, making evident (possibly approximate) symmetries that might be present. 

■ A testimony of the importance of this subject has been the number of works that have been devoted to it. From 
' the early works of York and Smarr 's', which proposed coordinate conditions through elliptic equations, to more 

■ recent works which present alternatives to choosing coordinates that could aid in the numerical simulation (see for 
I instance 0,|1,ISEjOi0|) considerable efforts have been invested towards defining useful coordinate conditions. In 

general, these conditions seek to minimize suitably defined quantities with the hope that these will, in turn, have a 
positive impact in the behavior of numerically evolved quantities. 
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I The proposed coordinate conditions are usually given in algebraic terms or through elliptic equations. The latter is 



expected when conditions requiring stationarity of variables- whose evolution is determined by hyperbolic equations- 
are imposed. When attempting to use such conditions, one faces the problem of dealing with a hyperbolic-elliptic 
^ ' system of equations which require appropriate boundary conditions. Until recently, convenient boundary conditions 
' for the main variables in spacetimes involving black holes were not sufficiently understood for generic situations when 
singularity excision was used [i^ • The complications associated with properly defining the elliptic side of the problem 
d ' at the analytical level coupled to the additional "extra cost" that solving these equations during the evolution has 
spurred a number of efforts aiming to circumvent both these issues. The idea has been to promote the elliptic equations 
to hyperbolic ones through "driver equations" These conditions aim to sidestep the cost issue and the need to 
impose boundary conditions at possible excision surfaces. 

Conditions based on this strategy have been employed in a number of works yielding much improved evolutions, 
most notably in BSSN codes where specific coordinate conditions are obtained by requiring the time variation of the 
(trace of the) connection variables be driven to zero (the so-called F-driver). The hyperbolicity analysis of the BSSN 
system with the F — driver conditions augmented by suitable advection terms has been presented in |36j | . It is shown 
that the system is strongly hyperbolic in this case, though unfortunately these augmented coordinate conditions do 
not necessarily freeze the F variables. It is then unclear wether these augmented conditions will have the same impact 
as the original ones in simulations and also if they share similar hyperbolicity properties. Therefore, there are a 
number of questions that remain open, namely (i) What is the true impact of adopting these 'driver' conditions on 
the hyperbolic properties of the complete (main variables plus gauge) system? (ii) Furthermore, are the conditions 
obtained sufficiently flexible to guarantee desirable properties such as to yield a convenient characteristic structure at 
boundaries? (iii) How must one extend the knowledge gained in the "F-freezing conditions" so as these can be thought 
as truly geometric expressions not tied to particular variables in the system -and hence useful to other formulations-? 
(iv) What is the freedom in the implementation of these conditions and their behavior in actual applications? 

In the present work, we examine these questions both in the theoretical and practical senses in order to draw 
conclusions applicable to most metric based formulations of Einstein equations by considering coordinate conditions 
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motivated from possible geometrical constructions. In particular, we concentrate on conditions defined either at a 
given hypersurface (and its embedding on the four-dimensional spacetime), or at the four-dimensional spacetime level. 
The former results into a set of elliptic equations which contains the well known minimal distortion/strain condition 
for the shift vector and the maximal conditions for the lapse while the latter gives hyperbolic equations related, in a 
sense, to the harmonic coordinates. 

With these conditions, we analyze the properties of the whole system of equations (coordinate conditions plus 
Einstein equations) where in the case of implementing the elliptic equations we promote them to hyperbolic ones 
via the "driver" approach. Additionally, we investigate possible difhculties that can be encountered when employing 
these coordinate conditions in conjunction with an excision strategy. 

Our analysis mainly concerns 3-1-1 metric formulations. That is, those based on the intrinsic metric and extrinsic 
curvature of spacelikc hypersurfaces defined by a foliation of the spacetime. The equations governing the future 
evolution of these variables are derived from the Einstein equations in the spacetime of interest and arc augmented 
by additional variables and check whether that the resulting system, coupled to the coordinate conditions, is at least 
strongly hyperbolic. 

As we will see, in all cases one obtains a system with a characteristic structure such that its eigenvectors couple 
gauge/coordinate variables to the main variables. This has strong implications for the system, since: 

• In the cases where singularity excision is to be used, the characteristics of the system must be such that they 
are completely outflow towards the excision boundary. This condition is fulfilled when, roughly /3" > a (with 
Z?" = /3*ni the projection of the shift along the spacelike unit normal to the excision surface Ui and a the lapse 
function). Since now the coordinate conditions are dynamical, extra care must be taken to monitor that it is 
fulfilled. 

• Since the characteristic modes of the system now mix coordinate and main variables, if the condition above is 
not satisfied, it is extremely difficult to provide boundary conditions to the gauge functions consistently. This is 
to be contrasted with the case where the elliptic conditions themselves are employed. Here boundary conditions 
for the gauge variables could be imposed so as to guarantee the outflow requirement is satisfied. 

Unfortunately, as we will describe in section III, commonly used conditions fail to satisfy some important desirable 
conditions which might have strong implications in numerical implementations. We point out how this can be 
avoided and the cost associated in doing so. Additionally, we analyze alternative conditions derived from considering 
approximate symmetries in the spacetime. This condition, called "harmonic almost-Killing equation" coupled to 
Einstein equations gives rise to a well behaved system where coordinates respond to (approximate) symmetries. We 
present simulations to investigate their usefulness within standard numerical relativity testbeds. 

In order to carry out several analysis presented in this work, it is necessary to adopt a given formulation of 
the equations. To this end we consider the strongly hyperbolic formulation presented in |l6j and the so-called Z4 
formulation. The former can be regarded as an 'augmented' ADM formulation with the addition of first order variables 
keeping the metric's gradient and suitable combination of constraints to the right hand sides. The latter is basically a 
covariant extension of the Einstein field equations, obtained by introducing a new four vector which is defined by 
its evolution equations. This way, the symmetrized covariant derivatives of this four vector are added to the Einstein 
Equations, that is 

i?M- + v^z, + v,z^ = 8 7r(r^, -ir^^,). (1) 

The solutions of the Einstein's solutions are recovered from the extended set when happens to be a Killing vector, 
that is 

V^Z, + V,Z^ . (2) 

For a generic spacetime, this happens just in the trivial case = 0, so true Einstein's solutions can be easily 
recognized. 



II. "IDEAL COORDINATES" 

Before presenting different possibilities for adopting coordinates we comment on what are arguably useful properties 
they should satisfy. There exist several discussions on what constitute requirements for good coordinate conditions 
(see, for instance, 1, 0|); these are based on the intuitive picture that useful coordinates, from the point of view 
of a numerical implementation, should: 
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• be free of artificial (coordinate) singularities; 

• take advantage of existing symmetries in the problem, whether approximate or exact. In particular if the 
spacetime is stationary, coordinate conditions should give rise to metric components explicitly time independent; 

• in the absence of symmetries, they should minimize the rate of change of either the metric or other appropriately 
defined geometrical quantities; 

• be 3-covariantly defined if possible (i.e. independent of coordinate changes within a given hypersurface) . 

The points above are important at a fundamental level in that they are hopefully satisfied irrespective of the 
formulation used or the particular problem under study. To these, further requirements might be added that refer 
more to specific applications like coordinates having: 

• suitable behavior near singularities. For instance, conditions yielding a convenient slicing of the spacetime. This 
could range from those that avoid the singularities altogether (like singularity avoiding conditions) to those that 
penetrate the possible horizon (in the case where excision techniques are to be applied). In the latter case, 
it important that the resulting characteristic structure be such that all variables are outflowing towards the 
excision region. 

• appropriate asymptotic behavior so that extraction of physically relevant quantities is facilitated and/or, related 
coordinate speeds are bounded so as to not have to deal with superluminal cases. 

Finally, the conditions adopted must be such that within the formulation of Einstein equations employed the well 
posedness of the underlying problem to be treated is guaranteed as the choice of coordinates is not decoupled from the 
issue of defining a well posed problem. For instance, the ADM formulation with analytically prescribed coordinate 
(lapse and shift) conditions is weakly hyperbolic (hence yielding an ill-posed problem) while with harmonic coordinates 
is symmetric hyperbolic. For a more general discussion of hyperbolic formulations sec for instance 

As mentioned, one of the main motivations when choosing coordinate conditions is that they should not introduce 
spurious "dynamics" in the evolution of the system. This motivation has lead throughout the years to the introduction 
of different conditions. When attempting to define such conditions an obvious difficulty is the need to do so in an 
three-covariant way so as to decouple coordinate effects to the true physical behavior of the spacetime. A way to do 
so was introduced by Smarr and York Q by constructing scalar quantities from the intrinsic and extrinsic curvatures 
of the spacetime and minimizing their variation with respect to the lapse function and shift vector. 

For instance, the 'strain' scalar defined in analogy with fluid dynamics, is defined as 

= = 2 (V(,/3j) - aK,j) (3) 

which can be used to construct a positive definite Lagrangian 

L = E., S*^' . (4) 

This geometrical object is used to construct a (non-negative) action which measures the distortion or strain of a given 
hypersurface. By minimizing this action with respect to the shift 

^ = 0, jL^d'x. (5) 

an elliptic equation (called "minimal strain") is obtained for the shift. When this equation is fulfilled, the rate of 
change of a suitable norm of the spatial metric will be minimized from one hypersurface to the next one. 

The minimal strain equation can be generalized by considering the action of the densitized metric 7^ 7^^ (with 
7 — det(7ij)), obtaining this way what we will call the "minimal densitized strain". This equation can be written 
simply as 

Vfe[S'^'' - A 7'^''' trE] = VklV'p'' + (3' - 2 A7'='V,„/3™ -2 a {K'^ - X-f^HrK)] = (6) 

so the minimal strain condition reduces to the choice A = 0. Another interesting case, called "minimal distortion" (in 
analogy with a related notion of elasticity) is recovered for the choice A = 1/3. 

These (elliptic) equations seem natural conditions for the shift in Numerical Relativity applications, since they 
satisfy the fundamental properties mentioned earlier. Notice that one could have also opted to minimize the action 
with respect to the lapse a. This yields an algebraic condition for it 0,13, though it is ill-defined for time-symmetric 
cases. An alternative strategy is to consider minimizing an scalar defined by KijK^^ with respect to the lapse. This 
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provides a fourth-order elliptic equation for it . We will refrain from considering it here as it will introduce further 
complications either at the computational cost level (if implementing the elliptic equation) or at a conceptual level 
when promoting the equation to a hyperbolic condition. 

As mentioned, when considering these conditions within an initial boundary value problem, assessing the well 
posedness becomes more involved as the system becomes of elliptic-hyperbolic nature. This, in particular, means the 
solution will strongly depend on the boundary data specified. Until recently, the lack of a well defined strategy to 
specify inner boundary conditions (say in the case where black hole excision is adopted), coupled to the extra cost 
associated with solving elliptic equations numerical induced considerable activity towards employing related conditions 
within a purely hyperbolic problem j43|. For these reasons, it has become customary to define associated hyperbolic 
equations to implement related conditions. Parabolic conditions could also be considered, but they certainly would 
not simplify the cost issue (as their Courant-Friedrich-Levy condition scales quadratically with the grid-spacing, and 
can be further regarded as an inefficient method to solve the elliptic equation itself via relaxation techniques). We 
will thus concentrate on hyperbolic conditions and analyze the implications they carry. 



III. HYPERBOLIC COORDINATE CONDITIONS 



We will restrict our analysis to a large family of hyperbolic coordinate conditions. Some of these are derived by 
simple relations that have been employed in the past while others are motivated by the minimization of geometrical 
quantities as described above. In this latter case, we follow recent works in that we will approach the problem 

here by adopting hyperbolic driver conditions to implement the equations. However, as opposed to these works, we 
will require that the coordinate conditions analyzed do indeed minimize the desired elliptic equations. That is, we 
will neither assume that considering other related equations having the same principal part but differing in the lower 
order terms will yield similarly behaved coordinates nor that suitable lower order terms can be added so as to obtain 
first order conditions. Although these assumptions can be, and have been, adopted in previous works, the conditions 
thus obtained are not guaranteed to satisfy the original sought-after geometrically motivated conditions. 

To be specific, we will consider conditions that can be written as 

dta = -o? Q , (7) 
dtP' = -aQ' (8) 

for suitably defined {Q, Q^} which we regard as the gauge conditions. These will be given by either algebraic or 
differential equations relating the Q-quantities with the other variables of the system, trying to fulfill as many of the 
desired requirements described in the previous section as possible. 

In what follows we consider separately three distinct cases which refer to the way the fields Q, are defined. We 
distinguish cases with the name 'algebraic' when Q or are directly defined, 'differential' when they obey evolution 
equations and 'semi-algebraic' in cases where an algebraic for one and a differential for the other is considered. 



A. Algebraic gauge conditions 

One of the simplest choices is an algebraic relation between the Q-quantities {Q, Q^} and the main variables of 
the system. In this definition are included the general gauge conditions proposed recently in [Toj and the subfamilies 
discussed in many other works (e.g., [1,0|)- The prototype of algebraic gauge conditions is given by the harmonic 
coordinates, which were introduced half a century ago to ensure the well posedness of the Einstein Equations [lll26l|. 
This is obvious as in this case, the principal part of Einstein equations for all components reduce to 

5"'(5,^),^, = Z.o. (9) 

Although the well posedness of the Cauchy problem is ensured this way for the evolution system, this coordinate 
choice does not fulfill "a priori" many of the properties of an ideal gauge condition in the absence of suitably defined 
sources or lower order terms. In particular, the freezing of the metric components in (almost) stationary spacetimes is 
by no means guaranteed. Another delicate issue is that the shift so-defined is not a three- vector and so need not reflect 
the symmetries in the problem. This can be seen more clearly by translating the harmonic coordinates condition to 
the 3-1-1 decomposition language: 

Q = -— dkh-ia + trK (10) 
a 

Q' = -— dkP' - a 7^=^ (5^7,fc - 9fe m V7 - dklna) . 
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The will not transform as a vector (except under linear transformations), so neither nor will be vectorial 
quantities during the evolution. 



B. Semialgebraic gauge conditions 

We will refer to as semialgebraic gauge conditions those that, keeping an algebraic relation for Q, allow for a 
differential definition of the . This way there is enough freedom to fix the shift with an exact geometric condition. 
Naturally, one could have done the opposite, i.e. an algebraic relation to while a differential one for Q Since the 
former is what is most commonly used in current applications and elliptic conditions for the lapse, that minimize the 
rate of change of the extrinsic curvature, is fourth-order we will concentrate on algebraic/differential conditions for 
the lapse/shift. 

For the lapse we propose a generalization of the harmonic coordinates which includes the Bona-Masso lapse condition 
[27| and its slight modification presented in 0, that can be written as: 

Q = -a—dklna + fia) (trK -2Q) (11) 

a 

where Q is added when considering the Z4 formulation of Einstein equations, otherwise this term must be dropped. 

The generalization consists on having added the parameter a to the above equations which determines wether the 
advection terms are included (a = 1) or not (a = 0). The "lapse speed" (ie, the speed of the eigenvectors associated to 
the lapse) will be fixed by this parameter in combination with the free function J [a). Notice also that the subfamily 
a — 1 reduces to several well studied cases depending on the expression for / = /(a): / = is the geodesic slicing, 
/ = 1 is the time harmonic slicing, / = 2/a is the "l+log" and so on. Additionally, the subfamily a = has been 
used successfully in the evolution of single BH 0. 

The rationale behind this generalization is that, as we will see later, the characteristic structure of the coordinate 
conditions adopted will have delicate, profound, differences depending on the values these parameters take. In addition 
to this generalization for the lapse condition, we consider a similar one to the shift condition which we will define by a 
suitable hyperbolic driver which seeks to satisfy the minimal distortion condition. We next revise how this condition 
is to be defined and further generalize it to include related options. 



Approximate geometric shift 

An elliptic condition can be imposed in a dynamical way through a parabolic or hyperbolic "evolution" equation. 
The former can be regarded as a standard relaxation way to obtain the solution of the elliptic equation while the 
latter drives the solution towards the desired one, in analogy with a damped oscillator. The first approach was used 
in |l5l | to convert the minimal distortion elliptic equations into time-dependent parabolic equations by means of the 
Hamilton- Jacobi method, that is, 

dtp' = aVk[^''' -^^^UrY.]. (12) 

The parameter a is characteristic of all the drivers, and determines the dissipation strength employed so that the 
solutions of the elliptic and the parabolic equations agree. For small values of ct, the shift is expected to tend slowly 
to the elliptic solution. At this point it is worth noticing that in the fully dynamical case, it is not completely clear 
that equation H12I) is actually a driver. The procedure is inspired in simple elliptic equations, like for instance the 
Laplace equation = 0. In this simple case the Hamilton- Jacobi method indeed provides a driver to the elliptic 
equation. In the case of Einstein Equations however, this conclusion is not immediate as the equations are highly 
coupled. In a "frozen variables" approximation, where all main variables are regarded as fixed, the driver condition 
does give rise to a solution satisfying the elliptic equations (up to suitable boundary conditions). In general, however, 
assessing this behavior is considerably more delicate. 

Nevertheless, current simulations indicate -at least for the cases considered- that the hyperbolic-driver conditions 
do give rise to reasonably well behaved solutions |^ as judged by monitoring the approximate fulfillment of the original 
elliptic condition that motivated the driver condition. These simulations implement a driver in such a way that some 
of the variables of the BSSN formalisms ,2^,22/1, the F*, are frozen at late times in black hole evolutions. Although 
they give rise to great improvement in the evolution of single black holes an head-on collisions, it is not clear whether 
they are successful due to the fact that they minimize particularly dehcate variables in the system ^ or due to their 
"proximity" (in a loose sense) to a minimal distortion condition. If the latter is the reason, it would indicate that this 
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condition could benefit other formulations. Unfortunately, to our knowledge, these conditions have been employed in 
practical applications only in the BSSN-based codes. 

In order to investigate the usefulness of the geometrically motivated condition we consider it within the driver 
approach generalized in the following way (Q3 equation): 

{dt -b£0)Q^ + Wk[ga {E>^' - A 7'=* trE)] ^ -a Q\ (13) 

Let us discuss in detail the differences between the standard gamma driver condition, as used in and l|13|l . First, a 
Lie term has been included in the Q3 equation with a free constant b. Although this Lie term does not come naturally 
from the "driver", we will see later that it affects the shift speed (that is, the speed of the eigenvectors associated to 
the shift) and it will be required in order to fulfill other requirements. The physical meaning, when b — I, would be 
that the driver is not along the time lines but along the normal lines to the space-like hypersurfaces. 

The second difference is that we employ a covariant derivative (with respect to the intrinsic metric of the hypersur- 
face) in (|13|l which is dictated by the minimal distortion condition. This additionally ensures the tensorial character 
of the equation, and so both the and the are now vectors, with the corresponding advantages. Finally, the 
parameter A has been kept in order to generalize the condition and adapt it to other formalisms that do not use the 
conformal decomposition. This way, it one can choose which densitized strain is going to be minimized during the 
evolution. 



1. Characteristic Analysis 

In order to analyze the structure of the system with the coordinate conditions adopted we must choose a particular 
formulation. Here we employ the Z4 formalism though we have checked that similar issues arise when employing the 
above coordinate conditions in the formulation presented in [T^ . 

The characteristic analysis of the gauge conditions I13f) with the Z4 formalism described in the next section 
shows that there are three clearly separated "gauge cones". We refer to them in this way to stress that these come 
about due to the coordinate conditions considered. However, their corresponding eigenvectors span not just the part 
of the Hilbert space corresponding to the lapse, shift and derivatives of this last one. Indeed, they have components 
both on the coordinates and main variables sectors which will have delicate consequences as we shall see later. These 
gauge cones can be grouped into three distinct entities: 

• Lapse cone, which propagates with speed — ^^Pn i \// + (^Y^)^ Pn ■ 



• Transversal shift cone, which propagates with speed —^^Pn ± y 5 ct^ + (^^)^ Pn 



• Longitudinal shift cone, which propagates with speed —^^Pn ± y 2 g (1 — A) a'^ + (^^)^ Pn- 

An analysis of the associated eigenvectors both for the Z4 and the formalism described in reveals that the full 
evolution system is strongly hyperbolic only if all the gauge speeds are different one from each other and different 
from the speed of light. Otherwise, there is a collapse of some of the eigenvalues and there is not a complete basis of 
eigenvectors, leading to a weakly hyperbolic system. 

We thus see the need to introduce the Lie terms (controlled with the parameters {a, b}) in equations ifTTlIT!^ which 
will provide sufficient fiexibility to obtain a well behaved system. For simplicity, let us focus on the condition for the 
lapse — the same discussion is also valid for the other gauge cones — . If the Lie term is included (a = 1), as the driver 
acts to minimize the dynamics along the normal line, the associated speeds are —P ± 1// a. This kind of structure 
allows for inflow coordinates where Pn > y/J a and all the lapse eigenvectors have negative speed. As mentioned, this 
requirement is crucial near a black hole horizon, where a standard practice is to excise the singularity by introducing 
an excision boundary. Here it is not known which boundary conditions to define and even how could be implemented 
if known ^4:^ . Another problem of this approach is the existence of so-called "sonic points" (in analogy with fluid 
dynamics) where the speed is zero (/3„ = \/J a). At these points there is a collapse of some of the gauge eigenvectors 
with some standing modes, and the system is weakly hyperbolic in one direction. 

On the other side, if the Lie term is not included (a = 0), the lapse speeds are — 1/3„ ± f a'^ + (%')^- In this case 
the driver is along the time lines and the evolution system is always strongly hyperbolic; the speeds/eigenspeeeds are 
such that only in the case a = there could be a collapse of eigenvectors. However, some of these eigenspeeds are 
such that, at an excision surface, will always describe incoming modes (i.e. towards the computational domain). This 
means that boundary conditions are to be specified for these modes somehow. Unfortunately, as these modes couple 
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coordinate and main evolution variables, one has to worry about how to provide suitable boundary conditions to the 
gauge functions and carry the evolution of the main variables without providing boundary conditions to them. This is 
a delicate problem in itself. Intuitively, one expects that boundary conditions are only required for the gauge functions; 
however, some of the main variables themselves depend on the gauge functions also (the extrinsic curvature). Hence, 
the issue of separating the gauge dependent component of the incoming modes must be clarified before proceeding 
this way. To do so requires considering constraint preserving boundary conditions which might be further complicated 
by the fact that the characteristic structure need not be constant along the excision surface. 

Summarizing, there is a tension between trying to obtain a minimizing prescription and ensure both strong hyper- 
bolicity of the system and that any suitably defined excision boundary is of outflow type. In hindsight it could be 
argued that this is a consequence of having tried to 'get away' without solving an elliptic equation -which does require 
boundary conditions at all boundaries- and solely deal with a hyperbolic equation where no boundary is required at 
the excision surface. 

It is then clear that the options are: (i) to stay at the elliptic (or related parabolic) level for the coordinate 
conditions; (ii) give up the symmetry seeking approach through driver conditions (at least in the problematic regions 
by suitably modifying the equations or by adding convenient lower order terms to the equations |3.ll | ) or (iii) consider 
a new set of options that aim to resolve the conflicts. 



C. Almost-Stationary Motions: the Q4 

An appealing alternative is to consider conditions derived by minimizing some suitably defined spacetime scalars. 
As it has been recently pointed out in the harmonic almost-Killing equation (HAKE) 

V^[^(^^'')-^(V.C).gn = (14) 

is a generalization of the Killing equation ^('^''') = whose solution space includes also the afline Killing vectors 
and, of course, its subfamily the nomothetic Killing vectors. For this reason, the covariant conservation law 114|l can 
provide a precise definition of the concept of approximate Killing vectors as solutions of the HAKE equation 114|l . 
This equation can be obtained from the standard variational principle (|SJ) with a Lagrangian L given by 

L = kp;^)&''^~\{'^-0^ • (15) 

Since the Lagrangian is non positive it is not possible to guarantee that extremizing the action will provide a solution 
that minimizes it. However, by suitably adding damping terms, the hope is that will indeed be the case. In such a 
case, the HAKE equation can be of great utility as a coordinate condition, because it is not only well adapted to the 
stationary spacetimes (a Killing vectors is a solution) but also "minimizes" the deviation from the stationary regime. 
In spacetimes with some (quasi) symmetry, it is expected that the congruence of time lines of the observers will be 
aligned during the evolution with the time (almost) Killing vector, avoiding this way spurious time dependence due 
to an unfortunate choice of coordinates. 

The physical meaning can be better understood by considering the adapted coordinates ^ — dt, where now = 
-C^Spiy — dtg^ij. The 4D Lagrangian (fT5|) can be written as 

L = S^, E^"^ - ig^, (S^5 5^*'). (16) 

which can be reinterpreted as a four-dimensional generalization of the positive definite 3D lagrangian Q). Following 
the analogy, the HAKE equation l|14() can be seen as the 4D generalization of the minimal densitized strain ((HJ. The 
main difference is that, since the HAKE equation considers also the time component of the spacetime, the structure 
of the resulting system is not elliptic anymore but hyperbolic. 

In the Z4 context there are Z-terms that must be included in the HAKE (|14|l in order to get a well posed problem. 
With these terms, the conservation law (|14(l can be written in different ways, like for instance 

V.[^/:c(V5 5n] = 2g^"^/:4Z. , (17) 

or, in adapted coordinates. 



(18) 
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Equation H17|) shows explicitly the tensorial character of the gauge condition, while the equation H18|) points out its 
relation with the harmonic coordinates, i.e. g'^P r^'^-p — 0. This gauge condition, which will be called Q4, is the 
closest to fulfill all the requirements; not only the shift but also the lapse is well adapted to stationary spacetimes, 
and if there is only an approximated symmetry, the coordinates are expected to adapt in order to minimize the rate of 
change of the metric. An additional property is that as a result of their construction, the gauge conditions obtained 
are also defined in a covariant way. 

The ambiguity of including or not the Lie terms in the Q-equations, introduced in the Q3 gauge, is not present 
here, where there is no choice: the Lie terms are actually there. As a consequence, the lack of strong hyperbolicity at 
sonic points appears again in the gauge cones, as it will be shown in the next section. 



IV. THE EVOLUTION SYSTEM: Z4 FORMALISM + Q4 GAUGE 

In order to study the hyperbolicity of the gauge conditions we have to consider them within the context of a specific 
formalism in order to get a closed set of equations that will constitute the evolution system. For concreteness we 
adopt the Z4 formalism, but similar results can be obtained with other formulations. 

Here the Z4 formalism and the Q4 gauge will be written down as an evolution system of (second order in space 
and first order in time) equations by means of the 3+1 decomposition. The characteristic structure of a fully first 
order version of this evolution system will be analyzed in detail, as well as how to pass from a second order system 
(in space) to a first order one without altering the structure of the eigenvectors. 

A. The Formalism : the (first order) Z4 system 

The four-dimensional equations can be written, by using the 3+1 decomposition, in the equivalent form j34j : 



{dt-Cp)-f,, - ~2aK,, (19) 

{dt - Cp) K.,j = -V^aj + a [R^J + V.Z^ + VjZ, - 2 + (tri^ - 2 9) K^j - + ^ {trS ~ r) 7.y] (20) 

[dt -Cf3)e = I + 2 VfeZ'= + {trK - 2 6) trK - tr{K^) - 2 Z'^Uk/a - 2 r] (21) 

{dt - Cp) Z, = a [Vj (K,^ ~ 5^ trK) + 9,9-2 K,' Z, - 9 - 5,] . (22) 

In order to convert the equations H19I22|I into a fully first order system, the spatial derivatives of the lapse, the shift 
and the intrinsic metric must be introduced as new independent quantities, that is 

A, = aana, = dkP\ Dk^^ = \^kl^J (23) 



and substituted everywhere. The evolution equations for these additional quantities can be computed easily taking 
the time derivative of the definition H23|l and permuting the time and spatial derivatives. Due to the commutativity 
of second spatial derivatives, we can add without any change in the solution space the constraints Cki = dkAi — diAk, 
Cik^ = diBk" — dkBi' and CkUj = diD^ij — dkDuj with free parameters {ca Cb Cd} to the evolution equations of the 
lapse, shift and intrinsic metric respectively. It the evolution equations for the metric components are defined in a 
general way as 

dta = ~a^ Q , dtP' = -aQ' , = ~2 a Q (24) 

then the evolution of their spatial derivatives, with the addition of the ordering constraints, would be 

dtA, + d,[a Q]~ca 13' Ch -0 (25) 
dtBk' + dk[a Q'] - Cb 13' Cik' - (26) 
^tDk^J + dk[a Q,j] - /?' Clk^J = (27) 

Here there is a delicate point; if one wants to preserve the same eigenvectors when passing from the second order 
system (in space) to the first order one, the choice {ca Cb Cd} = {111} is compulsory. Since we are interested on the 
physical solutions, which should not depend on the order of the (spatial derivatives of the) equations, this will be our 
choice from now on. 
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With this choice, a first order version of the Z4 formahsm can be written as a system of balance laws: 

(28) 
(29) 

Dk^J (30) 
(31) 
(32) 
(33) 

\\, = D\, - ^ (1 + (A/ + D,,'^) + U\ [A, + D, - (1 - E, - 2 Z,] (34) 

+ ij^- [A, + - (1 - e) ^» - 2 Z,] , 
being I^i = Dik' and i?i = D^^i. The non-zero source terms can be found in the Appendix. 



where 



dtA, - 


h di[ 


-/3' A, + 


(a g + /?" Am)] = 5/ - Bi Ai 


dtBk' - 


h di[ 




f (a + /3™ i?™')] - Bi' Bk' - Bi' 


dtDkij - 






+ 5\ [a Q„- + r D„,^j)] - Bfc' - 


dtKij - 


h dk\ 




+ a\\, ] = S{K,,) 


dtQ - 




[-/3'^ e + 


a {D^ -E^ -Z'')] = S{Q) 


dtZ, - 






- a {^K\ + 5\{trK - 6)} ] = S{Z,) 



B. The Q4 gauge 



The gauge condition can be written as a set of evolution equations for some gauge quantities by means of the 3+1 
decomposition, using either the covariant conservation law (|17|l or the non- vectorial "standing" equation p8|) . Of 
course, these equations (and their 3+1 forms) are completely equivalent, so one could be recovered from the other 
without any problem at the second order (in spatial derivatives) level. At the first order level this transition is not 
always so transparent when the ordering constraints C^i, C/fc' and CkUj are included; this is the reason to start from 
the appropriate version of the HAKE equation from the very beginning, at the four-dimensional level, to write then 
the most convenient 3+1 gauge equations. 

The 3+1 form of the conservation law (|17|l provides directly evolution equations for the tensor quantities {Q, Q*}, 
and it can be useful to take advantage of the symmetries of the problem. For instance, in spherical coordinates the 
vector is in general Q"^ — (Q*", , Q'^)- If the problem is also spherically symmetric, then only Q'^' and /3'' would 
have a non-trivial evolution equation, as opposed to what happens either with the harmonic coordinates (|10|l or the 
non- vectorial standing version (|18|l . Notice that the semialgebraic Q3 has also this vectorial character. 

On the other side, the 3+1 form of the "standing" version (|18|l reduces directly to evolution equations for some 
combinations of variables which follow an ODE, so they are directly standing modes of the system. This version is 
more convenient in general cases without symmetries in order to write the system in fully first order. The resulting 
standing modes, combinations of the Q-quantities with other variables of the system, hold always, so these eigenvectors 
are the same in both second and first order versions. The "standing" version (|18|) can be written, by means of the 
3+1 decomposition, as 

dtP = dt[a (g -trK + 2Q)+ (3^ Aj] = -2 [Q'^ ~ Q) ~ 2 a {A^ + Z^) (35) 

dtP' = dt [a Q' + (3^ + (2 E' ~ D' ^ A' + 2 Z')] = (36) 
2 a (a if/ - B/) + 2 (g^'= - 7-''^ Q) T'^k + 4 (Q'^ - 7'^ Q) Z, - a Q' [a {Q - trK) + A^] 

where the standing P-quantities have been defined. From these equations it is easy to see that the principal part is 
just the time derivative of the harmonic conditions (|10|) . 

As it was discussed previously, equations H35I36|I admit many different solutions. A convenient way to enforce 
the precise desired solution without unfavorably affecting the characteristic structure of an hyperbolic system was 
introduced in |35l |. The method consists in adding a source damping term that damps the solution to the desired one. 
In our case it would be: 

dtQ = ...-a {Q-T]trQ) (37) 
dtQ' - ...-<jQ' (38) 

where the dots stand now for all the original terms. Since only one vector can be constructed just by contracting 
the tensor with the normal lines, there can not be any ambiguity on the damping term for the equation of Q*. 
However, the two different scalars Q and trQ can be constructed from S^^, so all the combinations are included in 
the damping terms in (|37|) . Two special cases arise here: 
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• The first one would correspond to the choice 77 = 0; all the Q-quantities are driven to zero, so that one tries 
to minimize the rate of change of all the metric components. It is the default case, most suitable in physical 
situations in which we expect a stationary regime to be reached asymptotically. 

• The second case would correspond to the choice 77 = 1 ; the lapse equation H35|) is driven to the solution Q = trQ 
instead of Q = 0. This way, although the shift equation is still used for minimizing the intrinsic metric, the 
lapse just tries to follow the singularity avoidant condition dt{a/^) = 0. This can be useful in all cases in 
which singularity avoidance is required. Note that a stronger singularity avoidance behavior is expected when 
77 > 1. 



C. Characteristic structure 



The evolution system have the following 54 independent variables 

u = {a, P\ 7^,, K,j, e, Zfc, A,, B,\ Dkrj, P, P'} (39) 

where the {Q,Q^} can be written in the equations as function of the {P, P*} and other variables. The system is 
strongly hyperbolic if all the eigenvalues are real with a complete base of eigenvectors for any arbitrary direction 
and the symmetrizer can be shown to be smooth. We have not looked into this, though the analysis would follow the 
lines of those presented in 36, .S^l where the condition has been shown to hold. 

The analysis of the eingenvalue/eigenvector structure is more clear when the quantities (and the modes) are de- 
composed by projecting them into this specific direction. That way, for instance, a vector T!j would be separated into 
its longitudinal part T!„ = Tk and its transversal components Ta = Ti — Tn Ua- From now on we will use the 
indices {a,b,c,d} for the transverse components and n for the projection along n''. Using this notation, the list of 
the gauge-independent eigenvectors can be written as: 

• Standing modes: there are 10 eigenfields corresponding to the metric components with speed v — 

[«] , m , hi] ■ (40) 

• Normal modes: there 20 transversal components of the first order derivatives, propagating with speed v — — /?": 

[D,,,] , [AJ , [B,'] . (41) 



• Transverse light cone: there are 6 new independent eigenvectors propagating with light speed v = —f3n ± a that 
allow to recover {Kat, Dnab}, that is, 

Pa6<^' = [Kab- T^iBab + Pfca) ] ± [ A^ah - -^^^ {Dabn + Dfcan) ] • (42) 

2 a 2 



• Mixed light cone: there are 4 new independent eigenvectors propagating with light speed v = —(in ± a that 
allow to recover {Kna, Zg}, that is, 

ina'^' - [Kna] ± ^ I + i^ac" + (^ - 1) ^'ca - C D ann - 2 Z,, ] . (43) 



• Energy cone: there are 2 new independent eigenvectors propagating with light speed v = —(3n i a that allow to 
recover {Q, Zn}, that is, 

E^^^ = [e--B,^]±[Dnc'-D-,n-Zn] . (44) 

a 



Several comments are in order here. The eigenvectors are simple due to the choice {ca Cb c^} = {11 1}. Other 
cases, like the trivial one {c^ Cf, c^} = {0 0}, are considerably more complicated and do not have a complete basis 
of eigenvectors at the sonic points. Notice also that, up to here, the results are independent on the gauge condition. 
The remainder of the characteristic structure is dictated by the choice of coordinate conditions; in our case it is given 
by 
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• Lapse sector: The standing eigenvector [P], plus the cone spanned by the 2 new independent eigenvectors 
propagating with hght speed v = —fin ± a that allow to recover {i^„„, that is, 

G(±) = [P] + (« T [ {trK - 2 8) ± A„ ] (45) 

• Transversal shift sector: The standing eigenvector [Pa], plus the cones spanned by the 4 new independent 
eigenvectors propagating with light speed v = —(3n ± ot that allow to recover {Bna, Dnna}, that is, 

Sa^^'^ = [Pa] + {aTPn)[a{Aa+Da-2Ea-2 Za) ± [B^a + Ban) ] (46) 

• Longitudinal shift sector: The standing eigenvector [P„], plus the cone spanned by the 2 new independent 
eigenvectors propagating with light speed u = — /?„ ± a that allow to recover {i3„„, that is, 

Sn^^^ - [Pn] -{aTPn)[a Dn± (« tvK - trB) ] . (47) 

Note that at the sonic points (| /?„ |= a) one of the sign choices in the former equations coincides with one of 
the standing eigenfields {P, Pi}. Then, there is not a complete basis of eigenvectors and the system is just weakly 
hyperbolic there, as it happens with the Q3 gauge. However, these sonic points can be found only in the tachyonic 
coordinates regions, where 0^ > . Moreover, there will be missing eigenvectors there only for the particular 
directions in which | 0^nk |= ol. This is considerably less severe than the failure to achieve strong hyperbolicity 
for generic directions in the full computational domain. This so nic p oints issue arises also in the hydrodynamical 
equations and is often dealt with by a small amount of dissipation 38] . In the next section we will check numerically 
the behavior of the system. 



V. NUMERICAL RESULTS 

The previous gauge conditions will be tested in different periodic spacetimes which have been suggested |39l| as 
standard test-beds for Numerical Relativity codes. The numerical algorithm used is the standard method of lines |4Cl| | 
with centered second order discretizations of spatial derivatives and third order Runge-Kutta to evolve in time. We 
will focus on three different tests; the 'robust stability' test, in order to check the well posedness of the formalism. 
The gauge waves, in order to see the effect of the different gauge conditions in spacetimes with a time-like Killing 
vector. And the Gowdy waves, where there is no such time-like Killing vector. Usually we will compare the results of 
the (first order) Z4 formalism either with the harmonic coordinates (the evolution system will be called Z4harm) or 
the Q4 condition (the evolution system will be called ZQ4). The results for the zero shift case, which are identical in 
few cases to the Z4harm as we will show later, were already presented in [4]| . 



A. Robust stability 



Let us consider a small perturbation of Minkowski space-time which is generated by providing random initial data 
for every dynamical field in the system. The level of the random noise must be small enough to make sure that, as 
long as the implementation is stable, fields will remain at the linear regime even for a hundred crossing times (the time 
that a light ray will take to cross the longest way along the numerical domain) . This test is designed to experimentally 
assess the hyperbolicity of the evolution system by exciting high frequency modes and observing the overall behavior 
of the solution. As higher frequencies are allowed in the problem, for a strongly/symmetric hyperbolic systems the 
solution should be well behaved while this is not the case with weakly hyperbolic systems. 

The results of this test around the standard flat space-time /3' = are already well known from the analytical 
analysis for both the evolution systems Z4harm and ZQ4. Since both of them are strongly hyperbolic around — 0, 
all the norms remain constant during the simulation, decreasing slightly due to the inherent dissipation of the numerical 
scheme (no additional artificial dissipation has been added in the simulations). However, it can be useful to study 
numerically what happens at the sonic point = a, where the ZQ4 is strongly hyperbolic for all the directions but 
one, and check whether the ZQ4 system leads to a convergent solution. In order to study this scenario, we define both 
a and (5^ being 1 plus a small random perturbation. 

In order to see the expected behavior, we have plotted first in Fig. Qthe norm of trK of the Z4harm evolution 
system around the sonic point for three different space resolutions. Notice that, although we arc displaying one scalar 
quantity, the same behavior is observed by all the other norms. As it is expected in strongly hyperbolic systems, as 
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FIG. 1: Maximum norm of the trK for the harmonic gauge around the sonic point « a ~ 1 for three different resolutions. 
The slope of the norms remains constant independently of the resolution as expected on strongly hyperbolic systems. The 
simulation are performed in a cube of length L = 1 with 16, 32 and 64 points respectively. The time step is dt — 0.25 dx and 
no artificial dissipation has been added. 
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FIG. 2: The same plot that in Fig. Q but for the Q4 gauge. The slope of the norms, although they are growing in all the 
simulations, decreases as the resolution increases, approaching to the (constant) exact solution. This suggests that the system 
is well posed even at the sonic points. 

resolution is increased the numerical solution either should not grow or its growth should be lesser. Note also that 
the same kind of behavior is shown for both the Z4harm and the ZQ4 evolution systems around (3 = 0. 

A similar plot is presented in Fig. |21 for the ZQ4 evolution system, again for three different space resolutions. 
Although some of the variables (trK in the plot) show a growing norm, its slope decreases with resolution. So, 
although the profiles are different from the standard case shown in Fig. ^ the observed behavior is consistent wit a 
stable implementation, suggesting that the ZQ4 evolution system is not ill posed at the sonic points. 

B. The gauge waves 

We now consider again the Minkowski metric written in a non-trivial coordinates, obtained by performing a general 
conformal transformation to the t-x coordinates, that is. 



dx = 1/16 _ 




dx= 1/16 
dx = 1/32 
dx = 1/64 




20 40 60 80 100 



Propagation along the x axis can be simulated by considering a dependence like H(t, x) = h{x — t), so the exact time 
evolution (in these coordinates) will be just the "shifted" initial profile. We will use here a periodic smooth profile. 
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like a sine wave 



h{t, x) = 1 — A sin( 



2 TT 



(49) 



where d is the size of the x domain and A is the amphtude of the wave. AdditionaUy, we wiU take advantage of the 
periodicity of the initial profile to use periodic boundary conditions with d = I. 

In Fig.Othe norm of the strain Q^x is shown for both the Z4harm and the ZQ4 evolution systems. The evolution 
of the Z4harm is the same as the zero shift case described in ^] . In the ZQ4 case we have plotted two different 
cases, corresponding to different damping coefficients. The first one is with 77 = 1, so the time lines are driven to the 
condition Q = trQ. As it can be seen in the plot, the result is very similar to the harmonic case. The second case 
corresponds to the choice 77 = 0, where the time lines are driven to get aligned with one of the Minkowski time-like 
Killing vectors. The result shows the desired behavior; the observers behave in a way in which the metric components 
are explicit stationary, as it can be checked in Fig. ^ 

Different snapshots of the (non-trivial component of the) extrinsic curvature K^x, in Fig.|Sl shows that the evolution 
is almost frozen between 10 and 100 crossing times. Finally, a convergence test for the variable Qxx is performed in 
Fig. El The solution displays a decaying behavior, in which all resolutions match, until an asymptotic stage is reached 
(after around 30 crossing times), where the plots clearly converge to Qxx = 0. 
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FIG. 3: Norms of the non-trivial strain component Qxx for both the harmonic and the Q4 gauges with different values of the 
second damping parameter While the harmonic and the Q4 gauges with rj — 1 show a similar non-freezing behavior, the Q4 
with rj = actually minimizes the strain, driving the system to a stationary state. The initial amplitude of the gauge wave is 
^ = 0.1 and the simulations are performed in a channel of 50 x 5 x 5 points with length L = 1 in the longest direction. The 
time step is again dt = 0.25 dx and in this case some (small) Kreiss-Oliger artificial dissipation has been added in order to kill 
the high-frequency modes. 
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FIG. 4; The norms of the metric components for the ZQ4 evolution system with 7; = 0. After few crossing times all of them 
remain almost constant, implying a very small value of its time derivatives, the Q-quantities. 
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FIG. 5: The extrinsic curvature K^x in the x direction at different times for the same simulation that in Fig. 0| After 10 
crossing times there are not many changes in the profile. 

10 

1 

0.1 
0.01 
0.001 
1e-G4 
1e-05 
1e-G6 

10 20 30 40 50 60 70 80 90 100 

time 

FIG. 6: The non-trivial component Qxx is plotted for different resolutions [dx — 1/50, dx = 1/100, dx = 1/200) with the Z4Q 
with 1] = 0. All plots match during the transient decaying stage, until some minimum is reached. This minimum value can be 
seen to converge at a second order rate to Qxx = 0. 

C. The Gowdy waves 

Let us consider now the Gowdy spacetime, which describes a space-time containing plane polarized gravitational 
waves. The line element can be written as 

ds2 = gQ/2 („d^2 ^ j^2) ^ ^ ^^2 ^ ^-V ^y2) (gp) 

where the quantities Q and V are functions of t and z only and periodic in z, so that H5()(l can be implemented with 
periodic boundary conditions. Following 39], we will choose the particular case 

V = Jo{2nt) cos(27rz) (51) 
Q = 7rJo(27r)Ji(27r) - 27riJo(27rt)Ji(27r<)cos2(27rz) 
+ 2Tr^t^[Jo^(2TTt) + Ji2(27rt) - Jo^(27r) - Ji^{2tt)] 

(52) 

so that it is clear that the lapse function 

a = t-i/4 e2/4 (53) 



t=0 
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is constant everywhere at any time to at which Jo(27r<o) vanishes. In [S^l the initial shce t = to was chosen for the 
simulation of the collapse, where 27rto is the 20-th root of the Bessel function Jq, i.e. to — 9.88. 
Let us now perform the following time coordinate transformation 

t = ioe-^/^°, To = tjj/'*e2(*«)/4 ^ 472^ (54) 

so that the expanding line element l|50() is seen in the new time coordinate t as collapsing towards the t = singularity, 
which is approached only in the limit t — > 00. Notice that this singularity avoiding time coordinate r is not the proper 
time nor it does coincide with the number of crossing times, due to the collapse of the lapse. 

10 n 



QHr, = 0) 




time 



FIG. 7: Maximum norm of the lapse for both the harmonic and the Q4 gauges. After 100 crossing times the Q4 gauge with 
77 = gets too close to the singularity and crashes, while the other cases continue evolving until 1000 crossing times without 
problem. The simulations are performed in a channel of 5 x 5 x 50 points with length L = 1 in the longest direction. The time 
step is again dt — 0.25 dx and some amount of Kreiss-Oliger dissipation has been added. 
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FIG. 8: Maximum norm of the non-trivial component of the shift for the same simulation than in Fig. |7| The time scale has 
been changed from the linear to the logarithmic type in order to clarify the plot. Note that the Q4 gauge with rj — leads to 
a monotonic growing of the norm of the shift. In the other cases, we can see the strong oscillations of the harmonic case as 
contrasted with the smooth behavior of the Q4 gauge with rj = 1. 

As in the gauge waves test, the simulation is performed for the Z4harm and the ZQ4 evolution system with different 
values of rj. The maximum norm of the lapse a is plotted in Fig. [T] showing two different kinds of behavior 

• Singularity avoiding behavior. This is indicated by the collapse of the lapse, which can clearly be seen both in 
the harmonic gauge and in the ZQ4 case with rj = 1. This behavior is very similar to the zero shift case already 
described in [41|. We can see in Fig. |S1 that the rate of change of the shift is much slower in the ZQ4 case with 
77 = 1 than in the harmonic case, where strong time oscillations appear. This shows the freezing effect of the Q4 
gauge, when compared with harmonic simulations, even for singularity avoiding choices of the gauge damping 
parameters. 
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• Lapse freezing behavior. This is indicated by the absence of lapse collapse in the ZQ4 case with 77 = 0. Since 
the metric is collapsing, one gets close to the singularity in a finite amount of coordinate time and the code 
crashes. When translated in terms of proper time, however, all the simulations arrive approximately to the same 
point. We can see in Fig. [Sja sharp increase in (the norm of) the shift, which is trying to freeze the collapse by 
increasing the observers outgoing speed. 

It is worth to note here that this qualitative difference in the numerical simulations is triggered by the choice of the 
second damping parameter 77, without affecting the principal part of the original HAKE equation. 

The convergence of the solution for the ZQ4 system with 77 = 1 is shown in Fig.|^ where the Q scalar is plotted for 
three different resolutions. 
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FIG. 9: The G quantity is plotted for three different resolutions {dx — 1/50, dx = 1/100, dx = 1/200) with the Q4 gauge and 
= 1, showing a (second order) convergence to the exact zero solution. 



VI. CONCLUSIONS 

In this paper a large suite of hyperbolic gauge conditions has been studied with some detail, pointing out their 
advantages and possible problems. 

We have paid particular attention to conditions derived from geometrical scalars devised in order to minimize 
spurious coordinate effects. Our analysis reveal several important consequences applicable to these conditions and 
also all related ones (whose principal part coincide with those studied here, as the gamma driver condition): 

• Minimization of these quantities leads to a characteristic structure that yields inflow modes near the black hole 
excision surfaces. This implies that some kind of boundary condition is needed. However, as mentioned this is 
a delicate issue as main and gauge variables are intertwined in the inflow modes. 

• Related conditions, obtained by the addition of suitable advection terms to the equations, do resolve this issue 
but at a cost of bringing two more: the conditions do not necessarily minimize the sough-after scalars and there 
are surfaces where the system becomes weakly hyperbolic. 

Notice that there is a way to avoid most these problems altogether at the hyperbolic level by considering that a 
suitable first inte gral of the conditions does exist (which could be ensured by adding appropriate lower order terms 
to the equations) |lO|. However, the resulting conditions need not minimize the scalars and thus spurious coordinate 
effects might very well remain. 

As an alternative, a new coordinate condition, which has been introduced very recently, is used as a gauge pre- 
scription for Numerical Relativity applications. The main characteristic of this gauge condition (Q4) is that tries to 
"minimize" the deviation of the time lines from the time- like (quasi) Killing vectors, if there is one present on the 
space-time. The analogy with the 3D minimal strain condition is pointed out, and the evolution equations for the 
gauge quantities are written explicitly. The full list of eigenvectors is given, showing how to pass from a second order 
system (in space) to first order without changing the structure of the eigenvectors. In order to enforce the desired 
solution, some damping terms are included in the gauge equations, which allow for two kind of interesting alternatives. 
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The first one corresponds to freezing all metric components and it can be used when the space-time contains a (quasi) 
symmetry. The second one does not attempt to minimized the rate change of the lapse, but rather to drive it so that 
its rate of change is governed by the trace of the distortion. This provides a less restrictive alternative gauge condition 
for more general situations. 

Finally, some numerical experiments have been performed in order to check the properties of the Q4 gauge, the 
condition which appears as the most promising one within generic hyperbolic conditions derived in a geometrical 
way. First, with the robust stability test, it has been shown that the evolution system leads to solutions which are 
consistent with those of a well posed problem even at the sonic points, where the system is weakly hyperbolic just 
for some specific directions The gauge waves test is also employed to check the conditions, showing that the Q4 
gauge (for the choice t] = 0) indeed aligns the time lines with the time Killing vector, thus leading to a stationary 
state. The Gowdy waves test allows to further discriminate the effect of the damping parameters, leading to either 
a singularity avoidant or to a lapse freezing behavior (when the lapse is driven either by the condition Q — )■ trQ or 
Q — > 0, respectively). 
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Appendix: Sources of the Z4 evolution system 



S{K,j) - -K,, Bk" + K,k 5/ + Kjk + « {i (1 + e) h^fe T^, + ^{A, D, + A, A)] 

+ 1(1-0 [Ak D'^.,, - ^{A, (2 - A) + A, (2 E, - D,)} 

+ 2 {Dir™' D^rnj + Ar™ mi) " 2 {Dij'' + Dji'')] 

+ [Dk + Ak — 2 Zk) V^ij - T^mj r"'fei — {Ai Zj + Aj Zi) — 2 K^i Kkj 

+ (trK - 2 e) if,,} - 8 vr a [5,, - (-r + Sk'')] (A.l) 

SiZ,) = -Z, Bfc^ + Zk B,'' + a [A, {trK - 2 6) - Ak K\ - K\. V k^ + K\ {Dk - 2 Zk)] - 8 n a (A.2) 

5(9) = -e Bfc^ + ^[2Ak {D'^ -E''-2 Z'=) + Dk''' P".. - D''{Dk - 2 Zk) - K\ K^k 

+ trK {trK - 2 6)] - 8 tt a r 

(A.3) 
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